riwish::riwish(double innu, mat inS, icr_rands* inrands){
  //constructor

//uses Bartlett's decomposition to rapidly draw a wishart 
//and then inver it
  nu = innu;
//inS.print("inS:");
  S = inv(inS);
//S.print("S:");
  chol_S = trans(chol(S));
//chol_S.print("cholS");
//(chol_S*trans(chol_S)).print("chol_check:");
  rands = inrands;
  //cout<<"Nu:\t"<<nu<<endl;
 // gam_gens = new gamma_generator[nelem];
  for(int i=0; i < (int)S.n_rows; i++){
   //  gam_gens[i]=rands->rgamma_gen((nu-i+1)/beta);
     gam_gens.push_back(rands->rgamma_gen((nu-i+1)/2));
  }
  
}

////////////////////////////////////////////////////
mat riwish::draw(){
  double beta = 0.5; // for chisquare, beta =0.5
  mat output = zeros(S.n_rows, S.n_cols);
  for(int i=0; i < (int)S.n_rows; i++){
     //set the diagonals
    output(i,i) = sqrt(gam_gens[i]()/beta);
     //now the other bottom elements are all just normal deviates
     for(int j=(i+1); j<(int)S.n_cols;j++){
        output(j,i) = rands->randnorm();
     }
  }
  //take the inverse of L A AT LT
  return pinv(chol_S*output*output.t()*chol_S.t());
}

////////////////////////////////////////////////////
mat riwish::draw_noninv(){
  mat output = zeros(S.n_rows, S.n_cols);
  double beta = 0.5; // for chisquare, beta =0.5
  for(int i=0; i < (int)S.n_rows; i++){
     //set the diagonals
     //output(i,i) = sqrt( rands->rgamma_gen((nu-i+1)/beta)());
    output(i,i) = sqrt(gam_gens[i]()/beta);
     //now the other bottom elements are all just normal deviates
     for(int j=(i+1); j<(int)S.n_cols;j++){
        output(j,i) = rands->randnorm();
     }
  }
  return (chol_S*output*output.t()*chol_S.t());
}



////////////////////////////////////////////////////

riwish::~riwish(){
  //destructor
  gam_gens.clear();
}


